#include "core/fcc.hpp"

namespace FastMD
{
	
	static inline size_t fcc_meshes( Real x ) throw()
	{
		const size_t n = static_cast<size_t>( Anint(x) );
		return n <= 1 ? 1 : n ;
	}
	
	size_t FCC:: Count( U3D &cr, const size_t n, const V3D &box ) throw()
	{
		assert( box.x > 0 );
		assert( box.y > 0 );
		assert( box.z > 0 );
		assert( n > 0 );
		
		const Real vol1  = box.volume() / n;
		const Real alpha = __pow( __ONE / ( __FOUR * vol1 ), __ONE/__THREE );
		cr.x = fcc_meshes( alpha * box.x );
		cr.y = fcc_meshes( alpha * box.y );
		cr.z = fcc_meshes( alpha * box.z );
		return 4 * cr.items();
	}
	
	void  FCC:: Build( Vector &r, const U3D &cr, const V3D &box, const V3D &ibox, const Real noise_level ) throw()
	{
		const size_t n = r.vertices;
		assert( cr.x > 0 );
		assert( cr.y > 0 );
		assert( cr.z > 0 );
		assert( 4 * cr.items() == n );
		
		const Real dx    = box.x / cr.x; const Real dx2 = __HALF * dx;
		const Real dy    = box.y / cr.y; const Real dy2 = __HALF * dy;
		const Real dz    = box.z / cr.z; const Real dz2 = __HALF * dz;
		const V3D  noise(noise_level * dx, noise_level * dy , noise_level * dz);
		size_t     m = 0;
		/** SUBLATTICE A */
		
		r[0].x = 0;
		r[0].y = 0;
		r[0].z = 0;
		
		/** SUBLATTICE B **/
		
        r[1].x =  dx2;
        r[1].y =  dy2;
        r[1].z =  0;
		
		
		/** SUBLATTICE C **/
		
        r[2].x =  0;
        r[2].y =  dy2;
        r[2].z =  dz2;
		
		
		/** SUBLATTICE D **/
		
        r[3].x =  dx2;
		r[3].y =  0;
        r[3].z =  dz2;
		
		
		/** Repeat sublattices */
		for( size_t k=0; k < cr.z; ++k )
		{
			const Real z = k * dz;
			for( size_t j=0; j < cr.y; ++j )
			{
				const Real y = j * dy;
				for( size_t i=0; i < cr.x;  ++i )
				{
					const Real x = i * dx;
					for( size_t iref=0; iref < 4; ++ iref )
					{
						const size_t p = iref+m;
						assert( p < n );
						r[p].x = r[iref].x + x;
						r[p].y = r[iref].y + y;
						r[p].z = r[iref].z + z;
					}
					m += 4;
				}
			}
		}
		assert( m == n );
		
		/** create noise */
		V3D u;
		for( size_t i=0; i < n; ++i )
		{
			V3D &r_i = r[i];
			u.rand(1.0);
			r_i.x += u.x * noise.x;
			r_i.y += u.y * noise.y;
			r_i.z += u.z * noise.z;
			PBC_xyz(r_i, box, ibox);
		}
		
	}
	
}

